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ABSTRACT 

An  SQZ  algorithm  is  developed,  in  a  way  similar  to  that  of  Moler  and 
Stewart's  QZ  method,  for  handling  symmetric  A  and  B  where  B  is  an  ill- 
conditioned  positive  definite  matrix.   This  algorithm  preserves  symmetry, 
reduces  the  storage  requirements,  uses  less  time,  and  produces  only  real 
eigenvalues. 


Ill 
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CHAPTER  1 


INTRODUCTION 


In  many  applications,  such  as  those  in  physical  sciences,  the  solu- 
tion of  the  generalized  eigenvalue  problem  Ax  =  ABx  is  often  required, 
where  A  and  B  are  real  symmetric  matrices  and  B  is  positive  definite. 
There  exist  several  methods  for  solving  this  problem.   The  well-known 
Cholesky-Wilkinson  method  [l]  uses  Cholesky  factorization  of  B,  B  =  LL  , 
to  reduce  the  problem  into  standard  form.   However,  this  method  requires 
inverting  the  factors  of  B  which  may  lead  to  a  bad  solution  if  B  is  ill- 
conditioned.   For  a  nearly  singular  B,  Peters  and  Wilkinson  [k]    describe 
an  algorithm  which  approximates  the  null  space  of  B  and  removes  it  from  the 
problem  to  get  a  well-conditioned  problem.   This  method  involves  determining 
the  rank  of  B.   If  a  wrong  decision  is  made,  the  well-conditioned  eigen- 
values may  be  seriously  affected.   In  [6],  Fix  and  Heiberger  designed  an 
algorithm  which  is  a  variant  of  the  Peters-Wilkinson  method  [k]    for  nearly 
semidefinite  B,  i.e.  is  concerned  with  the  case  when  B  (or  A  and  B)  is  ill- 
conditioned  with  respect  to  inversion  and  B  is  permitted  to  be  positive 
semidefinite.   Peters  and  Wilkinson  also  describe  another  efficient  algo- 
rithm in  [2]  for  the  calculation  of  specified  eigenvalues  of  Ax  =  XBx  with 
band  symmetric  A  and  B,  the  latter  being  positive  definite.   In  their 
method,  every  eigenvalue  is  isolated  using  the  Sturm  sequence  property  of 
leading  principal  minors  of  A  -  XB  and  is  then  computed  accurately  using 
a  modified  version  of  successive  linear  interpolation.   Recently,  Moler 


and  Stewart  have  presented  the  QZ  method   [9]   which  was    designed  primarily 
for  nonsymmetric  matrices  A  and  B,    and  does  not   require   inversion  of  B. 
If  A  and  B  were   symmetric,   this   algorithm  destroys    symmetry  and  requires 
more   arithmetic   operations    and  storage  and  may  even  produce   complex  eigen- 
values. 

In  this  paper,   our  SQZ  algorithm  is   developed  in  a  similar  way  as 
the  QZ  for  handling  symmetric  A  and  B  where  B  is   an  ill-conditioned  posi- 
tive  definite  matrix.      This   algorithm  preserves   symmetry,    reduces  the 
storage  requirements,    uses  less  time,    and  produces   only  real   eigenvalues. 
Since  our  method  is   actually  a  symmetric   case  of  Moler  and  Stewart's  QZ 
method,   we  will   call  our  algorithm  "SQZ".      The  algorithm  is  based  on  the 
following  observations: 

1.  For  matrices  A  and  B   if  A  *   ZBZ     for  some  non-singular  matrix  Z, 
then  the  matrices   A  and  B  are   called  congruent.      If  A     =   ZA  Z 
and  B     =  ZB  Z     then  the   generalized  eigenvalue  problems 

A  x  =   AB  x  and  Ay  =   AB  y  have  the   same  eigenvalues,    and  the 
eigenvectors   are   related  by  x  =   Z  y. 

2.  For  matrices  A  and  B,    there   exist   unitary  matrices   U  and  V  such 

H  H 

that  both  A"    =  U  AV  and  B'   =  U  BV  are  upper  triangular.      The 

values   a'../b'..    are  the   eigenvalues    A.    in  Ax  =   ABx.       [3] 
n        n  l 

3.  If  A  and  B  are   symmetric  with  B  positive   definite,    then  there 
exists   a  matrix  U   satisfying  U  BU  =   I   such  that  A"   =  IjAU   is 
diagonal.      [3] 

The   algorithm  consists  of  the   following  four   stages: 
(i)   B  is   reduced  to   a  diagonal  matrix   (an   iterative  process),   while 
the  updated  A  is   still    symmetric.      This   stage  requires   only 


orthogonal  transformations. 

(ii)  A  is  reduced  to  the  tridiagonal  form  keeping  B  diagonal.   This 
stage  requires  both  orthogonal  and  elementary  transformations, 
(iii)  A  is  subjected  to  the  QR  transformations  with  elementary  trans- 
formations to  keep  B  diagonal,  an  iterative  process. 

(iv)  After  several  iterations  in  (iii),  A  approaches  the  diagonal 

form  and  the  eigenvalues  A.  will  be  given  by  a. ./b. .  if  b. .  ^   0. 

1  11   n     n 

If  a. .  ^  0  and  b..  =0,  then  we  will  have  an  infinite  eigenvalue 
n  n 

X..      If  both  a..  =  b..  =0,  then  any  scalar  can  be  an  eigenvalue 
i  n    n 

of  Ax  =  ABx. 
Stabilized  elementary  transformations  are  used  throughout  the  algo- 
rithm. 


CHAPTER    2 
SQZ  ALGORITHM 

2.1.      Theoretical  Basis 

Moler  and  Stewart's   QZ  algorithm  is   a  generalization  of  Francis'    QR 
method   [ll]  to   solve  the  problem  Ax  =   ABx  where  A  and  B  are   general   square 
matrices.      If  both  A  and  B  are  real   symmetric   matrices,   the  QZ  algorithm 
will   destroy   symmetry  and  hence   requires   more  time   and  storage  and  may 
produce  nonreal  eigenvalues.      This   is  not   economical   or  practical,    espe- 
cially  for  matrices   of  large   size.      As  we  mentioned  previously,    several 
algorithms   have  been   developed  for  dealing  with  the  above  problem.      How- 
ever,   it  was   mainly   assumed  that   B  in   addition  of  being  positive   definite 
is   also  well-conditioned.      In  this  paper  we  present   an  algorithm  to   deal 
economically  with  the  problem  when  B  is   ill-conditioned.      The   algorithm 
requires   much  less   storage  and  time,    and  produces   only  real  eigenvalues. 

Before  we  go  into  the  details  of  our  algorithm,  which  we  will  call 
SQZ,  let  us  present  the  idea  of  the  QR  method  for  the  eigenvalue  problem 
Cx  =   Ax, 

1.  Reduce  C  to   the  upper  Hessenberg  form  using  similarity 
transformations. 

2.  Find  an  origin  shift    A  using  the  roots   of  the  lower  right- 
hand  2x2  principal   submatrix  of  C. 

3.  Find  an  orthogonal  matrix  Q  such  that   Q(C   -   Al)   =  R,   where 
R   is   upper  triangular. 


k.      Let  C   =  QCQ   ,    then  the  matrix  C   is   upper  Hessenberg  again. 

5.  If  the  off-diagonal  elements   of  C   are  not  negligible,   then 
go  back  to  2. 

6.  The   eigenvalues  of  the  original  matrix  are  the   diagonal 
elements  of  C. 

If  the  original  matrix  C  is  real  symmetric,  it  is  first  reduced  to 
the  tridiagonal   form  which   is  preserved  by  the   QR  transformations. 

Moler  and  Stewart's  QZ  algorithm  is  motivated  by  the  QP  method  de- 
scribed above.      We  present  the   idea  of  this   algorithm: 

1.  Reduce  simultaneously  A  to  upper  Hessenberg  form  and  B  to 
triangular   form. 

2.  Find  the  origin  shift   using  the   roots   of  Ay  =   ABy,   where  A 
and  B  are  the  lower  2x2  principal   submat rices   of  A  and  B,    respectively. 

3.  Find  the  orthogonal  matrices   Q  and  Z,    such  that  QAZ   is   upper 
Hessenberg  and  Q(A  -   XB)    and  QBZ  are  both  upper  triangular  matrices. 

h.      Let   QAZ  be   denoted  by  A,    QBZ  be   denoted  by   B. 

5.  If  the  off-diagonal  elements   of  A  are  not  negligible,   then 

go  back  to   2. 

th  a'  ' 

6.  The   i        eigenvalue   is if  b..    ^   0. 

b.  .     n 
n 

The  idea  of  our  SQZ  is  similarly  presented  as: 

1.  Reduce  simultaneously  A  to  tridiagonal  matrix  and  B  to 
diagonal. 

2.  Find  the  origin  shift. 

3.  Find  a  transformation  L  such  that  LAL  is  tridiagonal, 
L(A  -  XB)  is  upper  triangular,  and  LBL  is  diagonal. 

h.      Let  LAL  and  LBL  be  denoted  by  A  and  B,  respectively. 


5.  If  the  off-diagonal  elements  of  A  are  not   negligible,   then 

go  back  to  2. 

th  a'  ' 

6.  The  i        eigenvalue   is  - —  if  b..    ^   0. 

ii 

To   simplify  the  explanation   in  our  SQZ   algorithm,   we   introduce  the 

following  notations: 

1.      Consider  the  real   symmetric  matrix, 


A  = 


T 


t*± 


*1 


!  0 

i   . 

t   i 

i   i 

i    i 

'    0 

t 

.th 

J        cc 

1 

0 


x2      0 0 


t 


.th        _ 
1        col, 


,th 

j        row 


.th 

l        row 


by  G,   we  mean  the   class  of  orthogonal   transformations   of  the   form, 


G.     .   = 


^  1 


.th 

i        row 


such  that  G.  .AG.  .  annihilates  elements  in  the  positions  (i,  j)  and  (j,  i), 
i.e.,  eliminates  x  in  the  matrix  A,  where  c  and  s  are  chosen  to  satisfy 


the  following: 


s   -c 


*1 


hence  we  have, 


*1 


r  '     r  ' 


nr   2 

•  V  X,  +  x,. 


r  =  sign(x1)«/x^ 
2.   Consider  the  matrix, 


D'  = 


'vX 


f        gr 


--x 


,th 

1   row 


hy  E,  we  mean  the  class  of  elementary  transformations  of  the  form, 


E.  = 

l 


»1    p 
0    1 


^1 


.th 

i   row 


;uch  that   E.D'E.    produces   0   in  the  positions    (i-1,    i)   and   (i,    i-l),    i.e., 


l        l 
annihilates  the  element   f  in  V    above  as    follows: 


'1       P 


0        1 


\1 


\      f 

f         K, 


^1  0 

P  1 


M 


.th 

l        row 


&,        o 


.th 

1        row 


-f 
where  p  =  — ,    and  g     =   g     +  pf . 

8p  -*-  -*- 


3.      By  S,   we  mean  the  class  of  diagonal  matrices   of  the   form, 


S.    = 

l 


.th 

l        row 


2.2.      SQZ  ALGORITHM 

We  will   start  "by   describing  the  algorithm  for  the   generalized  eigen- 
value problem 

Ax  =   XDx  (l) 

where  A   is    a  real  symmetric  matrix  and  D  =   diag(cL     ,    d      ,    ...,    d      )    in 

which   0   <   d, ,    <   d„   <    ...    <   d      .      Our  SQZ  algorithm  consists   of  two   steps, 
—     11   —     22  —  —     nn 

STEP   I — Reduce  A  to   the  tridiagonal    form,    keeping  D   diagonal: 
Let, 


A  = 


X         X X         x^  x? 

X  X XX  X 


*L 


X XX  X 


X_        X X  X 


and  D  = 


— l 


Consider  the  orthogonal  transformation, 


n,l 


s   -c 


From  (l),  we  have 

(G  ^G1   _)  (G   .X)  =   (G  .DG*  .  )  (G  _3fc) 
n,l  n,l    n,l        n,l  n,l    n,l 

or,        k'yi'    =  AD^x^ 

The  new  matrices  k'    and  D"   are  of  the  form, 


(2) 


A'    = 


x        x x        x        0 

X  X XXX 


X  X XXX 

0        x xxx 


,  v  = 


I 


Si   f 


'2  _ 


Now,   we   choose  the   elementary  transformation, 


E     = 
n 


1       P 
0        1 


—  f  t 

where  p  =  — ,  such  that  E  A^E  has  the  zeros  in  positions  (l,  n)  and  (n,  l) 
S2  n   n 

preserved,  and  E  D^E  is  a  diagonal  matrix.   Thus, 
n   n 

(E  G   iG1      E*)  (E-tG  ,*)  =  A(E  G  .DG1  ,  E* )  (E-tG   .5 ) 
n  n,l  n,l  n    n  n,l        n  n,l  n,l  n    n  n,l 


or,        k'SL'   =  AD'x' 
where  we  have, 


(3) 
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D'    = 


0 


0 


in  which 

d     -        nd 
n-l,n-l  nn        _ 

g     =   g     +  pf  =  *- >   0 

±  -1-  So 

In   general,    suppose  we  have  A  and  D  as    follows 


(h) 


A  = 


'i-1 


x 


h 

x^ 


.th 

J        col 


Xl      X2      ° 


X 


0 


t     J 


.th' 

i        col, 


.th 

J        row 


.th 

l        row 


D  = 


*1 


1 


.th 
-l        row 


where   T  is   a  tridiagonal  matrix  of  order   (j-l)   and  we  are   going  to 

annihilate  x      in  A,    so  we   apply  the  orthogonal  transformation  G.     .to    (l) 


and  obtain, 

(G.     .AG*    .)    (G.     .X)    =    (G.     .DG*     .)    (G.     .X) 

where  A"  =  G.  .AG.  .  and  D'  =  G.  .DG.  .  are  of  the  form, 


(5) 


11 


>1 


X  — 

I 

I 

0 


—  xT    0      0 0 


and 


g, 


^ 


respectively,   a   non-zero   element    f  is    introduced  in  positions    (i,    i-l) 

and    (i-l,    i)   of  D". 

Next,   we   apply  an   elementary  transformation  E.    to    (5)   and  obtain, 

(E.A'E*)    (eTV)    =      (E.D'E*)    (ETV)  (6) 

ill  ill 

where  E.    is   such  that   E.A^E.    has  the   same   zeros   as   in  A'   and  E.D^E.    is 
i  11  11 

diagonal.      The  elementary  transformations   E.    may  cause  numerical   instability 
if    |p|>l    (see  Wilkinson    [10],    p.    l6U ) .      For    |p|    =    I  —  |    to  be  less   than  or 


equal   to  1,   we  have, 


cs 


-1   < 


(Va2) 


<  1 


-     2,        2,      - 
s    d  +c    d_ 

dl  X2        s 

Let   r  =  —  >   0,    and  t   =  —  =  — ,   then 
d2  ^        c 

rt   +1 
Hence  the  region   in  the    (r,    |t|)-plane   in  which    |p|    >   1   is  the  shaded 
area  in  the   following  figure. 
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3  +  2/2. 


3  -   2^2 


v/2+1 


Suppose  that   x.^  ,    x0   and  d   ,    d     are   the   elements  which  yield    |p|    >   1   in 

A  and  D,    respectively.      In  this  case  we   apply  a  diagonal  transformation 

S.    in  which   sn    =  1   and  s,=  a  .      Therefore  S.AS.    and  S.DS.   yield, 
i  1  2  1111 


and 


r—              —* 

-'I 

1 

0 

*1 

Xl 

0 

a 

^2^ 

"2, 

1 

cT 

"dl 

0  " 

1 

— 

0 

di 

0 

0 

a 

0 

d2 

1 

0 

a 

0 

2. 

a  d 

thus,    t  =  at,    and  f  =  — — . 

a 

One  way  of   choosing  an  appropriate   scale   factor  a    is,    for  l6      <_| "b  |  <_  l6  , 

we  use  a  =  — ,    i.e.    t   =  1,   hence    |p|    <_  1.      For  the  other  values  of    |t|   we 
use  the  transformation, 


0 

/T 
1 
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hence  f  =  0  and  no  elementary  transformation  is  required.   Consider  the 

~  -2    i  i      2 

matrices  A  and  D,  where  l6   <  |t|  <  l6  , 


A  = 


"j-l 


x  —  x   x0  0  —  0 

i 


0 

0 

.th  II  ,    .tht  , 
J    col.   1    col. 


th 

j   row 


.th 

l   row 


D  = 


"2 


x  _J 


.th 

l    row 


and  the  corresponding  matrix,  (after  performing  orthogonal  transformation), 


D'  = 


52.. 


.th 

i   row 


in  which  —  >  1.   Consider  the  diagonal  transformation, 
So 


S.  = 

i 


Nl 


.th 
—  l   row 


1   Xl 
where  a  =  —  =  — ,  thus 
t   x2' 


l   l 


0 

|    x    — 

^ 

ax2 

0 

0 

0           |   x 

!ax2 

X 

!  o 

!    0 

lU 


and, 


l      1 


2. 

a   d 


2v. 


Therefore  the  multiplier  p  in   E.   will  he  less  than  or  equal   to  1   in   abso- 


lute  value, 


p  = 


f   |    i^V-V,    ,t2r-l 
_    l~l      p  2         '     '    2         '— 

s2  si+ci  t   r+1 


<  1» 


1  ,f    - 

since   r  >   0,    and  a   =  — .      For  those   cases  when     —     <   1,   we   choose   S.=I 

So  i 


f 


Now  if  we  define, 


j    j+2  j+1      n-1  n* 
where  V.  =  E.G.  .S.,  then  we  have, 

Z      -Z        ...Z0Z.Az!;Z*..Zt      Zt   =   T 

n-2  n-1  2  1     12  n-1  n 


and 


Z      _Z     n...ZQZnDZ^Z^...Zt   .Zt   =   D 
n-2  n-1  2  1     12  n-1  n 


j   =  1,    2,    ...,   n-2      (T) 

(tri diagonal) , 
(diagonal ) . 


That    is,   by  the  transformation   Z  =   Z        ...Z  Z    ,   we  will  he   dealing  with 
the  problem  Ty  =   ADy  instead  of  the  original   problem  Ax  =  XDx,   where 

T  =   ZAZ1 

D   =    ZDZ1 

X   =    Z    y. 

U  3 
To  reduce  A  to  T  and  D  to  D  we  require  roughly  —  n  multiplications  and 

12 

—  n      square   roots. 

STEP  II — Reduce  T  to  Diagonal   and  preserve  the   diagonal  D: 

This   is   an  iterative  process.      The  origin   shift    is  obtained  from  the 

-1/2     -1/2 
lowest  2x2  principal   submatrix  of  D  TD  .      Let  us   assume  that  the 


lowest   2x2  principal   submatrix  of  T  and  the   corresponding  2x2  principal 


15 


sub matrix  of  D  are  given  by 

B2   a2 


and 


ai    ° 


respectively.      The  shift   A    is   chosen   as   the   eigenvalue  of 


al        32 
B2       a2 


y 


di   ° 


y 


(8) 


or 


Tl        Y 


Y  T, 


y  =   Ay 


where  T      =   d     a 

T2  =   d"1a2 

Y  ..    ,-1/2-1/2 
-   d^        d2        32 

nl/2 

y  =  d      y 

which   is   closer  to  t  p.      Thus    A   is   given  by, 

x  ■  t2  +  r^ 


where 


q.  = 


*2  "  Tl 


=  ^77 


v  =  vq 

and  the  sign  in  q  +  v  is  chosen  to  yield, 

|  q  +_  v|  =  |  q  I  +  v. 
Wow,  let  the  matrices  in  Ty  =  ADy  be  given  by, 


(9) 


(10) 


(11) 
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T  = 


*1 
x„ 


2 

X  X 

XXX 

\  \  \ 

\  \  \ 

\  \  \ 

XXX 
X  X 


,   D  = 


X 

\ 


we   determine   an   orthogonal  transformation  G         to  produce 


s      -c 


x1   -   Xd1 


*1 

0 


(12) 


and  apply   it   to  T  and  D, 

(G2,1TC2,1>  (G2,ly)  =  UG2,iDG2,i)  (G2,iy) 

Thus   T    =  G       TG^        and  D'    =  G       DG^        are   of  the   form, 


T'    = 


X  X  + 


X  X 


+  X  X  X 


v  \    \ 

X  X 


\  N 


,    D'   = 


&1        x 

f  gr 


Next,  find  an  elementary  transformation  E  to  annihilate  f  but  does  not 
introduce  any  new  nonzero  elements  in  T" ,  thus  T  =  E  T^E  and  D  =  E  D'E 
are  of  the   form, 


T  = 


x       x1     x2 


XXX 


XX  X  X 

\        \         \ 
\  \  \ 


,    and  D 


IT 


Then  by  another  orthogonal  transformation  G    ,  we  annihilate  x  in  T  and 

introduce  again  a  nonzero  element  f  in  position  (3,  2)  and  (2,  3)  in  D, 

and  a  nonzero  element  in  positions  (U,  2)  and  (2,  h)    in  T.   By  another 

corresponding  elementary  transformation  E  ,  we  annihilate  the  element  f 

in  positions  (3,  2)  and  (2,  3)  in  D  again  and  do  not  introduce  any  new 

nonzero  elements  in  T.   That  is,  as  each  element  x~  in  positions  (i,  j) 

and  (j,  i)  is  annihilated  by  an  orthogonal  transformation  G.  . ,  it  pro- 

duces  a  nonzero  element  f  in  positions  (i,  i-l)  and  (i-1,  i)  in  D,  which 

is  immediately  annihilated  by  a  suitably  chosen  E..   In  the  same  way  as 

discussed  in  STEP  I,  a  corresponding  scaling  matrix  S.  may  be  required 

prior  to  performing  G.  .to  insure  that  the  future  multiplier  p  in  E.  will 

satisfy  |p|  <_  1  for  numerical  stability. 

The  process  continues  in  a  similar  way,  chasing  the  unwanted  nonzero 

elements  down  to  the  bottom,  right-hand  corner.   It  ends  with  M  =E  G    ^S 

n  n  n,n-2  n 

where  n  is  the  current  order  of  matrices  T  and  D,  and  we  have  obtained  a 
new  tridiagonal  and  a  diagonal  matrices  again.   If  this  process  is  applied 
iteratively  with  properly  chosen  origin  shifts,  there  result  sequences  of 
matrices  T^\  T(   ,  ...,  and  D(   ,  V        ,    ...,  satisfying 

T(i+D  =  z(±)^±)z(i)% ^  and  D(i+i)  =  ^lyiyi)^  (13) 

where      Z(i)  =  M(i)M(l> . .  .M,(i  V^V ±  \ 

n   n-1    4   3   2 

M{±)   =   EU)G?i}.  0SU),  for  j   =   3,  U,  ...,  n  (lU) 

and  the  special  transformation, 

(i)  _   (i)  (i)  (i) 
M2   -  E2  G2aS2 

which  is  produced  from  the  origin  shift.   The  matrix  T    will  approach 

diagonal  after  several  iterations,  replacing  the  negligible  subdiagonal 

elements  of  T    by  zeros.   So  the  eigenvalues  are  given  by, 
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a  . 

A      =  -J-,    if   3.    4   0 
J        3.  j 

where  T   "      =   diag[a.]    and  D   '       =   diag[  3  .  ] .      It   can  be  shown  that  one   iter- 

J  J 

ation  requires  roughly  22n  multiplications  (excluding  scaling)  and  n 
square  roots. 

To  replace  the  negligible  subdiagonal  elements  in  T  by  zeros,  we 
use  the  following  two  criteria,  consider  any  2*2  corresponding  diagonal 
submatrices  in  T  and  D,  assume  that  these  are  as  given  in  (8).   From  (10) 
we  see  that  &     is  negligible  if, 

l)  |y|  5.  eM  ( I  T-i  I  +l  To  I  )  is  satisfied,  where  e  =  machine  epsilon,  or 


2)  both 


Tl 


:  l6e„,  and 


xj  -    M 


<    e.,  are   satisfied. 


i  T   t     '    —     M 
12 

The  eigenvalues  of 

_Tl        Y    1 

y  =  *y 

J  T2_ 

are   given  by, 


2  2 

X      -    A(t     +    t_)    +    Tiljl    -  -* — )    =   0 
1  2  12  t    t2 

2 
thus,    if  the  term  -r— ^ r  <_  eM,    i.e.    negligible,   then  our  origin   shift  will 

be  just   the   same   as   if  6      =  0   in  T. 

One   special   case  may  be   disposed  of  immediately,    since  we   assume  that 
the  matrix  D   is   positive   definite   and  its   diagonal   elements   are   in  ascend- 
ing order,    so   it  may  happen  that  D  has   cL      so   small   that   dividing  by   it 
may   cause  overflow.      In  this   case  we  use  the   following  criterion:      If 

d..  _    <   e.,  is   satisfied,    then   replace   dn  _    by   0.      Thus   T  and  D   are  of  the   form, 
11  M  '  ^  11     J 
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T  = 


*1      X2 


XXX 

\  \  \ 

XXX 


,    and  D   = 


22 


M 


nn 

In  the  above  situation  we  use  an  elementary  transformation  E  to  anni- 
hilate the  element  x~  in  positions  (2,  l)  and  (l,  2)  in  T,  then  the 

t  t 

matrices  T'  =  E  TE  and  D'  =  E  DE  are  given  by, 


T'  = 


x 

0 

0 

X           X 

X          \      NX 

\    \ 

N  x     Nx 

and  D'   = 

x  J 


-x. 


t  2 

where  the  multiplier  of  En  is  p  =  . 

2         x^ 


If  |p|  >  1,  we  just  scale  both  T 


-m 


and  D  by  the  matrix  S  =  diag(q,  1,  ...,  l)  with  q  =  l6  ,  where  m  is  chosen 
as  the  smallest  integer  such  that 


*1{ 


<  1 


i°gjp 


i.e., 


-]  < 


log  l6   — 
e 


rn. 


Now  we  go  back  to  our  original  generalized  eigenvalue  problem  Ax  =  XBx. 
As  we  have  mentioned  in  section  2.1  above,  that  before  we  deal  with  the 
problem  Ax  =  ADX,  we  have  to  diagonalize  the  matrix  B  first.   This  diagonaliza- 
tion  can  be  done  by  applying  our  SQZ  algorithm  to  the  standard  eigenvalue 
problem  Bx  =  yx  to  solve  for  the  eigenvalues  and  eigenvectors  of  B.   Then 
those  eigenvalues  are  the  diagonal  elements  of  D.   Of  course,  every  trans- 
formation which  is  applied  to  B  in  solving  Bx  =  yx  must  also  be  applied 


to  the  matrix  A,  i.e.  if  ZBZ  =  D,  then  A  =  ZAZ  .   Therefore,  the  gen- 


20 

— "t      ~      — t 
eralized  eigenvalue  problem  Ax  =  ABx  is  reduced  to  A(Z   x)  =  AD(Z  x), 

i.e.  AX  =  ADX,  where  X  =  Z   x.   The  problem  Ax  =  ADX  can  then  be  solved 

as  we  have  discussed  above. 

Note  that  only  orthogonal  transformations  are  applied  in  the  pro- 
cess of  solving  for  the  eigenvalues  and  eigenvectors  of  B. 

It  is  interesting  to  compare  the  number  of  operations  required  to 
solve  the  eigenvalue  problem  Ax  =  ABx  by  the  SQZ  and  QZ  algorithms  where 
A  and  B  are  real  symmetric  and  B  is  positive  definite.   We  will  assume 
that  the  QZ  will  be  applied  without  modification  so  in  reducing  B  to  the 
upper  triangular  form  the  symmetry  of  A  is  destroyed.   For  the  SQZ  al- 
gorithm, assuming  that  we  need  both  the  eigenvalues  and  eigenvectors,  the 
operation  count  can  be  divided  as  follows: 

(i)   Bx  =  yx: 

h      3 

1.  In  reducing  B  to  tridiagonal  form  we  require  roughly  —  n 

1   2 
multiplications  and  -r- n  square  roots. 

2.  Reducing  the  resulting  tridiagonal  matrix  to  the  diagonal 
form  essentially  is  an  iterative  process  that  requires  17n 
multiplications  and  n  square  roots  per  iteration. 

3.  Forming  the  corresponding  matrix  A  consists  of  the  following 
three  parts: 

(a)  Saving  all  the  orthogonal  transformations  in  1  above 

3 
requires  roughly  2n  multiplications. 

(b)  Saving  all  the  orthogonal  transformations  in  2  above 

2 
requires  roughly  i+n  multiplications  per  iteration. 

(c)  Forming  XAX  =  A,  where  X  is  the  orthogonal  transforma- 

3 
tion  formed  in  (a)  and  (b)  above,  requires  roughly  2n 

multiplications . 
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(ii)     Ax  =   ADx: 

1.  In  reducing  A  to  the  tridiagonal   form  T  and  preserving  the 

Ij.     3  12 

diagonal  D,   we  require    —  n     multiplications   and  —  n      square 

roots. 

2.  Reducing  the  resulting  tridiagonal  matrix  T  to  the  diagonal 
form  and  preserving  the  diagonal  D  essentially  is  an  iterative 
process  again.   In  this  reduction  step  we  require  22n  multi- 
plications and  n  square  roots  per  iteration. 

3.  Forming  the  eigenvector  matrix  Z  consists  of  the  following 
two  parts: 

(a)  Saving  all  the  orthogonal  transformations  in  1  requires 

3 
roughly  2n  multiplications. 

(h )  Saving  all  the  orthogonal  transformations  in  2  requires 

2 

roughly  Un  multiplications  per  iteration. 

The  operation  counts  given  above  can  be  summarized  as  follows: 


TABLE  2-1 

OPERATION  COUNT  OF  REDUCTION  STEP 
(For  SQZ  it  consists  of  the  steps  (i) 
(ii)-l,  and  (ii)-3.a  given  above.) 


without  eigenvectors 

with  eigenvectors 

* 

r~ 

* 

r~ 

SQZ 
QZ 

20  3,,,  2, 

—  n  +  kn   k 

IT  3 
3~n 

2 
n  +  nk 

2 

n 

26  3   ,,  2. 
.—  n  +  4n  k 

^3  3 
6~ n 

2 
n  +  nk 

2 
n 

(*) 


(*) 


k  is  the  number  of  iterations  required  in  step  (i)-2. 
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TABLE  2-2 


OPERATION  COUNT  FOR   ONE   ITERATION 
(In  SQZ  we  have  Ty=ADy  where  T  is 
tridiagonal   and  D   diagonal, 
in   QZ  we  have  Hy=ARy  where  H  is   an 
upper  Hessenberg  and  R  upper  triangular. ) 


without  eigenvectors 

with  eigenvectors 

* 

r~ 

* 

r~ 

SQZ 

QZ 

22n 
13n2 

n 
3n 

kn2 

2 
21n 

n 
3n 

If  k  =   2n   and  assuming  that  we  need  two   iterations  per  eigenvalue   X   in 
Table   3-2,    then  SQZ  will   require  a  total   of: 


Without   Eigenvectors 
88     3 


multiplications  while  QZ  will   require: 

Without   Eigenvectors 

190     3 
—  n 


With   Eigenvectors 

1U8     3 
—  n 


With  Eigenvectors 
295  3 


multiplications,  where  we  assumed  that  n  is  so  large  that  n  and  n  are 

3 
negligible  compared  to  n  . 
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CHAPTER  3 
NUMERICAL  RESULTS 

The  SQZ  algorithm  has  been  implemented  in  a  FORTRAN  program.   The 
program  is  written  for  finding  the  solutions  of  the  generalized  eigenproblem 
Ax  =  ABx  for  real  symmetric  matrices.   Generally,  we  call  our  SQZ  program 
twice  for  the  full  symmetric  matrices  A  and  B.   The  first  call  of  SQZ  reduces 
B  to  a  diagonal  matrix  and  updates  the  symmetric  matrix  A.   If  B  is  already 
in  diagonal  form,  then  the  first  call  of  SQZ  program  is  not  necessary.   The 
second  call  of  SQZ  program  finds  the  eigenvalues  and,  if  required,  the  eigen- 
vectors of  the  system.   Between  the  first  and  the  second  calls,  we  rearrange 
the  diagonal  elements  of  B  such  that  they  are  in  ascending  order  to  yield 
better  results.   All  the  examples  presented  in  this  paper  have  been  run  on 
the  IBM  360/75  at  the  University  of  Illinois. 

In  all  the  examples  we  have  seen  so  far,  our  SQZ  algorithms  require 
roughly  1.3  iterations  per  eigenvalue.   Also,  we  tend  to  find  the  smaller 
eigenvalues  first  so  we  can  guarantee  that  the  smaller  eigenvalues  computed 
by  SQZ  will  be  as  accurate  as  possible. 

We  compare  some  results  obtained  by  SQZ,  QZ,  and  TRED2  and  IMTQL2  [7] 

-1/2  -1/2 
applied  to  B    AB    .   The  program  IMTQL2  or  program  IMTQL1  may  yield 

better  results  if  we  arrange  the  diagonal  elements  of  T,  the  tridiagonal 

-1/2  -1/2 
matrix  obtained  from  B    AB    ,  in  ascending  order.   Thus  for  those  ex- 
amples which  produce  the  diagonal  elements  of  T  in  descending  order,  we 
flip  the  matrix  T  before  we  call  IMTQL1  or  IMTQL2  to  guarantee  better  re- 


2k 


suits.   Also,  for  some  other  examples,  we  compare  some  results  obtained  by 
SQZ,  QZ,  and  Kahan  and  Varan's  program  "RECURSECTION"  [5]. 

In  the  SQZ  program  we  actually  combine  the  orthogonal  transformation 
G.  .  and  the  corresponding  elementary  transformation  E.  as  a  single  trans- 
formation  as  follows: 


E.G.  . 
1  i,J 


1   P 
0   1 


1 


\ 


s   -c 


1 


Ncl   si 
s     -c\ 


1 


where  cl  =  c  +  ps,  and  si  =  s  -  pc . 

Furthermore,  we  print  out  the  normalized  eigenvectors  in  our  program 

output.   And  the  relative  residual  is  computed  as, 

I  3.  Ax  -  ct.Bxl  I 
1  '  1      1   1  '00 


(1% 


+  la.  I   |b|   ) 


,th 


where  a .  is  the  i   diagonal  element  of  the  final  diagonal  form  of  the  ma- 
trix A  and  3.  is  the  i   diagonal  element  of  the  corresponding  diagonal  ma- 
th ai 
trix  B.   The  i   eigenvalue  is  given  by  — —  .   Numerous  examples  had  been 

3i 
tested  using  SQZ  and  most  of  them  gave  results  agreeing  to  at  least  lk   sig- 
nificant figures  with  QZ,  TRED2  and  IMTQL2,  or  RECURSECTION  [5]  (their  ALGOL 
program  has  been  translated  into  FORTRAN  on  the  IBM  360/75  at  the  University 
of  Illinois  by  the  author).   Tne  following  carefully  chosen  examples,  how- 
ever, indicate  some  interesting  disagreement  in  the  results  obtained  by  SQZ, 
and  QZ.   For  those  underlined  eigenvalues  we  find  that  there  is  more  agree- 
ment between  SQZ  and  TRED2  and  IMTQL2  (or  TRED1  and  IMTQL1 ) ,  while  QZ  gives 
a  poorer  result. 
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Example  1: 

2 

1 

3 

1+ 

A  = 

1 
3 

-3 
1 

1 
6 

5 
-2 

k 

5 

-2 

-1 

_0  [.  Q 

,    B  =  diag(lO      ,   10"    ,   1,   10    ) 


SQZ: 


(b)   QZ: 


Result 

2.00005  00362  50999 

-3.U9991  96525  38813 

1.571U2  5781+5  8771+8 

-5.16363  I+878I  31513 


-3.50000 

71 1+25 

3903 

-5.16363 

1+8717 

6210 

1.56603 

53712 

1+899 

(c)    TRED1 

-5.16363 

I+878I 

311+78 

and 
IMTQL1 : 

1.571U2 

5781+5 

87753 

■3.1+9991     96525     38812 
2.00005     00362     51001 


xi  0 
1 
xio 

xlOC 

xlO 


-7 


xlO 
xlO 
xlO( 


-7 


xlO 

xlO( 
1 
xlO 

xlO 


-7 


Number  of 

Residual 

Iterations 

10-21 

0 

10-18 

1 

10-16 

0 

ID"19 

2 

10-16 

0 

io-1T 

0 

10-12 

0 

ID"11 

2 

2 

1 

1 

0 

Example   2: 

"  6 

1+ 

1+ 

l 

A  = 

1+ 
1+ 

6 
1 

1 
6 

1+ 
1+ 

1 

1+ 

1+ 

6 

ft  li  ft 

,    B  =   diag(lO      ,   10"    ,   1,    10    ) 
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(a)   SQZ: 


(b)    QZ: 


(c)    TRED1 
and 
IMTQL1: 


Result 

6.00026  67081  U6831 

3.33326  85379  I+287O 

2.U9993  75718  79956 

-7.1+9999  96999  70010 


3.3331+1  66729  2031 

-7-^9999  96139  U273 

2.101+86  86976  6099 

-7-1+9999  96999  69996 

2.1+9993  73718  79958 

3.33326  85379  U2865 

6.00026  67081  1+6830 


*10 
1 
xlO 

xlO 

xlO 


0 


xlO 

xlO 
xlO( 


-8 


xlO 

xlO 

1 
xlO 

xlO 


-8 
0 


Number  of 

Residual 

Iterations 

ID-22 

0 

10-19 

0 

ID"17 

1 

lO-21 

2 

io-ig 

0 

lO"" 

0 

10-1° 

0 

10-1° 

2 

2 

1 

1 

0 

Example   3: 

2 

l.uitei 

A  = 

1.I+1U21 

-5.98802     -i+.  31621 
-i+.  31621     -2.51908 



1.38862 

Result 

(a)   SQZ: 

2.00005 

00358      89225 

-3.U9991 

161+17      62911+ 

1.571^3 

621+71+      91697 

-5.16353 

76268      1+9121+ 

(b)   QZ: 

2.00005 

00358        8923 

-3.1+9991 

16U17        6292 

-5.16353 

7621+3        11+1+8 

1.56990 

86I+36       5hki 

1.38862 
1+.  1+8198 

xlO8 


xlO 

xlO( 

xlO 


-7 


xlO 

} 
xlO 

xlO 

xlO( 


-7 


o 

B  =  diag(lO   ,  1.998206x10 
,-2 


-1+ 


9.338338x10      ,   I.67I+716XIO1) 


Number  of 

Residual 

Iterations 

10-25 

0 

10-18 

1 

10-15 

0 

10-18 

2 

io-3i 

0 

10-27 

0 

10-12 

0 

10-12 

3 
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(c)   IMTQL1:        -5.16353  76268  U9086  xio 

1.571^3  62^7^  91707  ><10( 

-3.U9991  161417  6291^  xlO 

2.00005  00358  89226  xlO^ 


-7 


Example  h' 


where 


A  = 


P  = 


ft  L.  —ft 

and  D  =   diag(8><10    ,    7,   10"    ,    2x10      ) 


2 

1 

3 

1» 

1 

-3 

1 

5 

,    B  = 

PtDP 

3 

1 

6 

-2 

wu 

5 

-2 

-1 

r~ 

1 
2 

1 
2 

1 
2 

1 
2 

1 
"2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
"2 

1 
2 

1 
_   2 

1 
2 

1 
~2 

1 
2 

(a)   SQZ: 


(b)   QZ: 


Result 

2.28911  60739  92863 

-2.26278  38980  361lU9 

-8.57121  08967  89169 

1. 3lU8l  U8112  63813 


xlO 
xlO 

xlO 
xlO 


-1 

-8 


-1.08900   87389    8385    xlO 


7 


-8.57121     085^7       2372 
1.31U81     U8112       6382 


-1 
-8 


Number  of 

Residual 

10-1^ 

Iterations 
0 

10-1^ 

1 

10-17 

1 

10-18 

1 

10-16 

0 

10-16 

0 

10-16 

0 

10-16 

3 
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(c)    TRED2 
and 
IMTQL2: 


1.31U81  U8112  63815  xio" 

■8.57121  OT213  7^638  xio 

8.57^69  90230  68577  xlO- 

•8.56361  98936  36599  xlO' 


-1 


Example   5: 


where 


A  = 


P  = 


1 

-3 

1 

5 

3 

1 

6 

-2 

1+ 

5 

-2 

-1 

1 

2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

1 
2 

,    B  =   P  DP, 


and  D  =  diag(2xl010,    7*10    ,    9X102,    3><10"10). 


(a)  SQZ: 


(b)  QZ: 


Result 


1.U8710   73189   3^508   xio' 
-1.U8699   ^9075   32303    xio' 


-8.57039   538U2   65270    xlO 

5.25869  12660  U5T56   xio 


-7 
-10 


-1.00251  66ohk      8551   xio' 

t 

1.00258   55293    UU70    xlO' 


■8.57039   53838    5853 
5.25869   12662    9592 


-7 
-10 


Number  of 

Residual 

It 

erations 

10-^ 

0 

io-lU 

1 

ID"15 

0 

10-16 

2 

10-16 

0 

10-16 

0 

10-16 

0 

ID"17 

2 

(c)    TRED2 
and 
IMTQL2: 
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5.25869      12662      95921  xlO 

-8.57039     53838     65502  xio 

9.M399     35712     6910U  xioJ 

-9.1+1+396     95971     92862  xio] 


-10 
-7 


For  these  underlined  eigenvalues,  there  is  only  one  significant  figure  agree- 
ment among  all  three  methods  while  the  other  eigenvalues  agree  as  expected. 


Example  6: 


A  = 


10 
1 


8L 


-9     1 

1   -10    1 


"181 

1   9 

1 


B  =  diag(l,  2,  3,  ...,  37,  e^  Eg,  £3,  e^) 


where               e  =  l6  x  10 

e  =  16  x  10 

e  =  16  x  10 

e,  =  16  x  10 


-10 
-11 
-12 
-13 


1 
10 


1+lxl+l 
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(a)  SQZ: 


6.25685 

99337 

201+76x10 

-3.1+5956 

99738 

86960x10 

5.56U10 

58UT8 

29690X1011 

-2.39999 

08392 

l+l+o68xio_1 

h. 93162 

79161+ 

70630xl010 

-2.91631 

661+91 

l+l+512xl0_1 

1+. 28820 

23318 

53898xl09 

-3.12U83 

683U9 

82l57xlO_1 

1.82853 

1 1+710 

09l+05xlO_1 

-2.66666 

06521 

50588xio_1 

1.1*5237 

26693 

97532xiO_1 

-2.1U285 

70906 

63091X10"1 

1.15VT9 

Uli+07 

U9982xi0_1 

-1.5381+6 

15381 

95596xio-1 

8.83UOT 

59855 

l6290xiO~2 

-8.33333 

33333 

25l+l+Oxio~2 

6.06101 

16889 

29363xl0"2 

-3.3221+3 

68998 

8923l*xl0_1 

3.12500 

61+277 

981+50xiO~2 

3.1+0U89 

59l+l6 

90256xlO_1 

5.36158 

2U036 

535UI4X10"1 

1.00000 

00000 

00102xl0_1 

3.33333 

3331+8 

390U5xlO~2 

2.22222 

22222 

2l+368xlO_1 

h. 85622 

5125^+ 

9U727X10"1 

3.75000 

00000 

1+5261X10"1 

6.89655 

172U1 

1+U9U0X10-2 

5.711+28 

57152 

66886xio_1 

I+.I+360U 

16099 

62082xlO_1 

8.33333 

33551+ 

70959xlO_1 

1.0T1U2 

85713 

30581X10"1 

1.20000 

00527 

1+8301x10° 

l+.ll+lll 

19392 

21581+xio"1 

1.75000 

13591+ 

80318x10° 

3.89^89 

1U106 

951+OUxio"1 

2.66670 

62771 

92957x10° 

I.I+81I+8 

1I+76I+ 

91+751xlO~1 

U.5011+3 

8221+1 

27761x10° 

1.92307 

66966 

o6680xio-1 

1.00898 

09703 

9U306X101 

3.5^661+ 

61165 

055l+2xio-1 

(b)  The  corresponding  results  in  QZ  and  RECUESECTION  agree  to  at  least  Ik 
significant  figures  with  the  results  from  SQZ  shown  above,  except  for 
the  two  smallest  underlined  eigenvalues  as  expected.   The  largest  resid- 
ual computed  in  SQZ  is  10~   while  that  of  the  QZ  is  10~   ,  the  differ- 
ence is  due  to  the  fact  that  QZ  uses  orthogonal  transformations  all  the 
way  while  SQZ  uses  also  elementary  transformations  which  may  affect  the 
accuracy  of  the  computed  eigenvectors. 
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APPENDIX: 
FORTRAN  PROGRAM 
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IMPLICIT  PEAL-fl    (A-H,C-!I 50JJ3J1 

LOGICAL  WANTX  S0Z0002 
DIMENSION    A  \* 1,41)  ,  *(< ,  1  ,41 ) ,X K 1 .Al  I ,  S  49  (  82  3  J ,  At  f  A(  A  1 J ,  BET  A  (  ^1) SQ7CO03 

I3IHF.NSIJN   A  DIM  I,  HUN  1 1,  l>  I -VI)  i  lit* 1411  .W1U1J S'.'/J.13v 

FOUIVALENCE     ( A ,B » , ( D, BETA  J  SQ70005 

cccccccccccccccccccccccccccccccccccccccccccccccrcccccccccccccccccccccccc   sozo^oe 

T rHlb    MKUUhA"    IS   WMlllbN    IU    S>iLVb     IHb    NbNbPALIZbU    b 1 UbNVALUb    lJ»JHLb", b'.'/JJiH 

C       A*X=LAMPPA*R*X,    WHERE    A    ANO    6    APE    REAL    SYMMETRIC    MATRICES    AND    B    IS    AM    SQZ0308 
C       ILL-CONDITIONED    POSITIVE    DEFINITE    MATRIX.    THIS    PROGRAM    IS    A    SYMMETPIC    SQ70009 

T CASL    Oh    MULbH    AND"    b  !  b'.JAM  I  « i>    HZ    PklH^AM    \J4±    MU!     II    P^b5>b<Vtb    SYW-'.b  !  k  Y  ,     i'JZmiO 

C       REDUCES    THE    STORAGE    REQUIREMENTS,    USES    LESS    TIME,     AND    PRODUCES    ONLY         S0Z0011 
C       REAL    EIGENVALUES.      SO70012 

_c 5>U7UD13 

C      USAGE    :  SQZOOIh 

C  REAL*8    A«  ND,ND»,^IND,ND)  ,X(ND,ND)  ,  S  AB  (  ND<M  NO- 1  ) /2 ), AD<ND) , POIND)     SQ70315 

T FbAL*B    ALFAINU)  ,-^blAINJJ,  dl NUI ,W1 I NU J  ,1  IcH  IUUJ STTTTFTr 

C  EQUIVALENCE     (  A  ,  B  )  ,  (  D,  BETA  )  SQZ0017 

C  .  S"?0318 

t : STTTWrr 

C  .  S070020 

C        CALL  SOHND,  N,B,D,  EPS,.  TRUE.  ,X,'£2  , 1  E»R  ,  £•♦,  ITER,  IFLAG)  SGZ0321 

-r . 507*022 

C                         .  SOZ0023 

C  .  S*Z0024 

T CALL    S0ZJKD,K,A,D,EPitWANTV.,X,E3,IEP,Bf  E/;,ITErMFLAP.) S1TZ7T72T" 

C                         .  SQZ0026 

C  . SCZ0327 

~C : ~~ SU7TT77B" 

C                     END  SQZ0029 

C  WHERE    THE    PARAMFTE°S    USED    APE    :   S^ZOP30 

T INPUT    : SOZflnjl 

C  A,    B  -    CONTAIN    THE    N    BY    N    REAL    SYMMETRIC    MATRICES    AND   USE  SOZ0032 

C  THE    SAME    ARRAY    STCRAGE     INDEPENDENCY    IN    l-ST    AMD  SOZ0033 

T 2-ND   CALLING  Gf    SU^SuUTINE    S07. $0  70034 

C  SAB  -    IS    A    VECTOR    KHICH    HAS    ND*(ND-l)/2    ELEMENTS    AND  SOZ0335 

C  IS    USED    TG    SAVE    THE    ORIGINAL    INFORMATION    CF    CNE    CF  SOZ0036 

T THE   TvO   SY^METAIC    KATPICE5    A   AMI   B   whILF  "HE   OTHER 5T70TTT 

C  IS    STORED    IN    THEIR    SHAPEING    ARRAY    STORAGE.  SGZ0038 

C  N  -    ORDER    OF    tHE    INPUT   MATRICES.  S070O39 

"C ND" MJ«bEp   OT  F.CwS   OR   Columns'    OP  AAA  Ays'    SPECIFIED. Sqzoc'6 

CD  -    CONTAINS    THE    DIAGCNiL    ELEMENTS    OF    IDENTITY    MATRIX    FCR  SO70041 

C l-ST    CALL    CF    SO*,     AND    THE    EIGENVALUES    OF    "ATPIX    6    FOR  SQZO'W 

~C"  2-ND    CALL    CF    SO*.    THIS    VECTOR    USPS    THE    SAME    ARRAY  SQZ00-»3 

C  STORAGE   CF    VECTO°    nETA(ND).  SQZ004* 

C AD -    CONTAINS    THE    OIAG°NiL     ELEMENTS     OF    ORIGINAL     DATA    OF     A.  SOZHQ',? 

~C"  BD  -    CONTAINS    THE    DIAGONAL    ELEMENTS     Or    ORIGINAL     DATA    CF    B.  SQ7T>t 

C  X  -    CONTAINS    N    RY    N    IDENTITY    MATPIX    FCR     l-ST    CALL    CF    SQZ,  SOZr>047 

_C AND    THF    TRANSFORATION    "ATRI^CS     IN    PI A GONAL I Z I MG       B  SQ7Q0~0 

C  FOF     2-ND    CALL    OF     SO?.  SQ7"Ot9 

C  EPS  -    IS    THE    RELATIVE    PRECISION    CF    ELEMENTS    OF       A    AND    fl.  SQ70050 

_C I T  L  A  G -    INPIC^E*    THE    r  F^U  ic  EMENT    OF    PERFORMING    ELFMENTj-CY  SOZOOrl 

C  TRANSFORMATIONS,  SGZQ152 

C  IFLAG  =  Q    -    INDICATES    ELEMENTARY    TR  ANSFCJR  VAT  I CNS    A*  E  SQZ0053 

_C pEQIII  CF0, ; . SOZOOf- 

C  1FLAG=1    -    INDICATES    NO    bLE"cNTA^Y    7*  ANSFOF*AT  I  CNS  SQ70n?3 

C  APE    FEOHIPED.  SOZOO"^ 

_C WANTX -    IT'DIC^TES    THE    c  EP'JI  c  E  *'F  NT    CF    COMPUTING    E  T  CEVV  EC  TOR  S  ,  Sn*03c7 

C  WANTX=.TKUb.    -    MEANS    EIGENVECTORS    A^E    REQUIRED,  SUZ0053 

C  WANTX=.FALSF.-    MFANS    EIGENVECTORS    A-;E   NOT    REQUIRED.  SQ7O059 

C  OUTPUT    ; SOZQ  HO 

C  ALFA, BETA    -    ARE    VECTORS    AND   CONTAIN    THh     INFORMATION    OF    EIGEN-  SOZ^QSl 


3U 


T VALUES.    THE    I-TH   EIGENVALUE    IS   GIVEN   flV , : SflZflOo2 

C  LAMRDA(I)    =    ALFACI )/«ETA(I I.  SQZ0063 

C  A,    B  -   THE    STRICTLY    UPPER    TP I  ANGLE  <  I .  E .     A(I,J)«0.    FOP     I>«J  )       SQZOOo-V 

T CUNTAIMS    THL   LI-klGlMAL    U^PKHAMuN    AHO    TZ   SAVED    TflFPtF    S«n3o* 

C  GHOUT    THE    COMPUTATION    3UT    THE    LOWER    TRIANGLE    IS  SQ7TOo6 

C  PVERWklTTEN    AND    ONLY    THE    DIAGONAL    WHICH    CONTAINS    THE       SQZO^o? 

T VALUED   DF    ALF-M11    IS   NEEDED   TIJ   iJUTFUT,    THE    57C  ICTLV S0731od 

C  LOWEP    TRIANGLE     IS    DESTROYED.  SQ7DC69 

CD  -    CONTAINS    THE    VALUES    or     RETA(I).  SQZOOZO 

T K CCMAIK5    I  hE    tIRENVcfT  rss    IP    REQUIRED. SnZfl.171 

C  ITEP  -    ITEPIII    CONTAINS    THE    NUMBER    CF    ITERATIONS    NEEDED    FOR       SOZOD72 

C  FINDING    THE     I-TH    EIGENVALUE.  SOZ0073 

T TFRU IS   TROUBLE    INDICATE, 537^174 

C  IERR»J    -     INDICATES    THE    J-TH    EIGENVALUE    HAS    NOT    BEEN    SQ70175 

C  DETERMINED    AFTE&    CEPTAIN    AMOUNT    CF    NUMBER       SOZ0O76 

T Uh    ITFOATIMMSi  Wt    St'    J  =  JU    IN    I  HIS    PF-C'CAMI    SOfiMJ 

C                  £J                      -    INDICATES    SPECIAL    RETURN    POINT,    J,     FOR    COMMENTS    ON  SOZ0078 

C  TROUBLES. SPZ0J79 

T  STT77TTJTT 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC    Snzo08i 

c  sozooa^ 

T IBM  S1ANUAUU  HXUP  bUbkUUIINb  LI-  EkkUW  420 !    I    <i21'6    : SH70^d3 

CALL  ERRSET<207, 256, -1,1, 0,208)  S0Z008^ 

WANTXs.TRUE. SO70O85 

EPS=2.22D-U SOZHOdo 

ND*M  ■■...„  SOZ0387 

1  CONTINUE  S0Z0088 

READ(b,10,bND=99tl    K ?Cn09 

WPITE(6,9)    N  SOZT09D 

9    FORM AT (1  HI,'    N    =',13//, '    P    =  ') SCZ0091 

rc-TTT^rnTD $0703-12 

DO    6    1=1, N  S0ZDT93 

READ(5,90O)  (B<  I,  J)  ,J  =  1,N)  SOZ009^ 

6    WPITF(6,10nO)    (a(I  ,JJ,J  =  1,NJ 507  M45 

900    FCP*AT(8C10.2)  SOZ0096 

AKCRM=0.  SQ70097 

BnnPM=o. snzoj^a 

IJ  =  0  S0Z0O99 

DO  11  1=1, N  S070100 

D(I  J  =  1.D0 : ?T7TTTr 

PNI=0.  SO7O102 

DO  20  J=l  ,N S0ZO1P3 

rni=rni*dabS1&I  i,j|  j  swtttw 

c initial  i 7f  x  :  sq701p3 

20  X(  I  ,  J)=n. SOZri06 

PNORM=DM,i*l  (RNr  P.M,?.N  I  )  SO  ^  31  M 

RD(I)=P«I,I)  S0Z0108 

X(  I  ,1  )=1.D0 SQ/OI09 

IP1  =  I*1  SQ'0110 

IF(I.CE.N)  C-P  TO  11  SQZOlll 

DO  21  J  =  IP1,N S0'P112 

1  J= I J  +  l  SOZ3113 

21  sap(  i.n  =  «(  I ,  J)  snzmu 

11   CONTINUE SQ711  15 

1000  FTPMATI'  ',^12011.3)  SQZ>H16 

NM1=N-1  S0ZP117 

N2  =  N/2 , S^'Ulfci 

IFL£r,  =  i  S0^119 

CALL     SOZ(ND,M,P,D,CDS,  .TRUE.,X  ,r.?  ,1  ERR  ,  J>,  I  TEf. ,  IFL  AG  )  SO7D120 

IFl  Af~  =  P S0Z0121 

PRINT    80P  SQZU122 
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800  FORMAT!/,'  A  > 
DO  24  1-1, N 
24    0(  I  l  =  B(  1 ,  1 1 


5AJ0123 
SQ!012'> 
S0Z1125 


1)1'    <!b    1*1, N 

5Q73125 

READ<5,900)       (A(  I,J),J»1,NI 

SOZ0127 

WRITE (6, 1000)!  A( I ,J) ,  J»1,N» 

SQ70128 

2b        Cl'NTlNUE 

5070129 

DO    27    1=1, N 

S0Z0130 

ANI»0. 

S0Z0131 

AUl ll>AlltI I 

1>U/<U132 

DO    26    J*1,N 

SO/0133 

26    AM  =  ANI+DABS(M  I  ,JJ  ) 

SQZ0134 

21    ANUKM  =  DMAX1  UNnKP.Z.NTJ 

SO7013F 

CALL    OPDEP!N,ND,D,X) 

SQZU36 

PRINT    700,     (D(  I)  ,1=1, N) 

SOZ1137 

/UO    RKMAI|/,»    EIGENVALUES    GF    UJTKIX      B '«•/ ,  (vX  ,  !PbU2S>.l6  H 

S071138 

C    PEPFCP*    TPANSFO"MATICNS    ON    MATRIX    A    : 

SOZ0139 

CALL    XTCX(ND,N, A,X, AD.W1) 

$070140 

CALL    SUZINU.Nti  ,P,  LPS,  WANT  X,X,  E3"  ,  I  EPR"  ,&*  ,  IT  EP  ,  I  FLAG  } 

SO  7  "J  141 

DO    35    1=1, N 

SQZ0W2 

35    ALf  M  I) =A< I ,  I) 

S0'0143 

NMl=N-I 

SP7J14A 

IF(  .NCT.WANTX)    GO    TO    32 

S0701<5 

C RELOAD      B      f.      SAVE      A    : 

SQZ01<-6 

1J  =  0 

snziHJ/ 

DO   40    I-1,N 

S0Z0148 

IP1=I*1 

SO/0149 

ihi.r.n  i.iv  ; l  *i 

SO71150 

DO   ''0    J  =  IP1,N 

SQZ0151 

IJ=IJ+1 

SO70152 

T  =  AU,J) 

SOZ0153 

BII, J)=SAB( IJ) 

$070154 

40    SAMIJ)  =  T 

S07.T155 

41    CONTINUE                               "" 

sozoist 

CALL    XTCX(ND,N ,B,X,BQ,W1) 

SOZ015  7 

32    CONTINUE 

SPZ0158 

DO    997    1=1, N 

S0Z0159 

IF     (BFTAU).EO.O.)    GO    TO    999 

S0Z0160 

EIG=ALFA(I  J  /PETM  I  ) 

S0ZD161 

WF  HE  16,993)    BIG, ALFA!  I  )  ,  «ETAI  I  ) 

SOZ0162 

998    FORMAT!/'    LAMBDA*    •  , 1 PD26  .16 , / '     • , • ALF A , BET  4  =  ' , 1P2D26 . 1 6 , 

SCZ0163 

♦           /3X, •EIGENVECTOR*' J 

S07016A 

GO    TO    150 

SOZ0165 

999    WFITF(6,99C  )    Al  F  A  (  I  )  ,  BETA! I  ) 

SOZ0166 

995     FrRM#T{/«     ','LAMqDa=    *  *     ^FINITE    »*'/,'     • , ' ALFA , 3ETA =• ,  1P2D26  .  16  , 

Sy/nL-7 

♦          /3X,  'f  ir,FNVFCTCP=«  ) 

SOZOlori 

150    CONTINUE 

S0Z0169 

IF( .NCT.WANTX)     GO    TO    997 

SO/117  0 

Wl ( ! )=0SC*7!DARS(B< 1,1))) 

SOZ0171 

DO    33    K=l,N 

SOZ0172 

33    XI  K  « I  )=X(K,I  )/Wl  (I  ) 

SO/T173 

PRINT       9<^,     (  XI  J,  I  I  ,  J  =  l  ,N> 

SO/.")  17' 

99A    FCPMAT( •»• ,<X, 1PD26.16) 

SPZ0175 

FN    =     0. 

SO/  -)17^ 

KJ  =  0 

SOZ0177 

DO    c5     K=1,N 

SOZ1178 

KMl=K-l 

S0711 79 

KP  1=K+1 

SOZ0180 

SA=AO(K)*XCK,I ) 

sozoiai 

S«=     FP(K)*X(K,I ) 

SO/0132 

ir  (k.fo.n)  r-r  t."  «i 


SO  Z -J  Id  3 
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dp  50  J=KP1,H 

KJ=KJ*1 

S£  =  SA*SABe<J»'X<  J,I  ) 


bH*i>B*K(  K  ,  J  |  ^  V  C  J  ,  I  J 

50  CONTINUE 

51  CONTINUE 


507 0ld4 
SO  70 135 

snzoiafc 

SO Z 018  7 
S0711J8 
SO/0189 

SOZ 3191 
S0Z0192 


kc  =  kj-2  .no*  ir-K) — ~ 

IF(K.EQ.l)    GO    TO    5  3 
DO    c2    JJ=1,K«1 


~J=Kri-JJ  +  l 

SA=SA*SAR(KC)*X( J,I ) 
KC=KC-( N-K+JJ) 


bl     bB=bB*bU,K  J'XI  J,i) 

53    CONTINUE 

P    =    OAPS(RETA(  I)»SA-ALFMI)*SB) 


TTmrrrr 

SQ/D194 
SO/0193 
SUZJNg 
S0Z01lJ7 
SO  70 198 
S'JZ'.'ISV 
SOZ0200 
SO/0231 
~W7"77T2~ 
SOZ0203 
SQ™20'» 

SO/02  03 
SO/0207 
SOZ3238 
S0ZO209 
SO/0210 
S0ZJ211 
S0*0212 
SO/3213 
5QM214 
S07021c 
SO/0216 
TOToTTT 


kn  =  i;vaxnhfhN) 

55  CONTINUE 

IF    (RN.NE.O.) 


1 KN    =    >-N/IUmS(^klf  I  1  )  J»ANlkM+lJAm;(ALFMlH-*!31MJaM) 

PPIMT    208, PN 
997    CONTINUE 


2~0B    H'M'ATI    3X,  "«b*>n)UAL»'  ,1PU26.16///) 
GO    TO    1 
2    PPINT    221  ,     IEFP 


2^1    f-UKHAT(////l3X,«««IHfc    ',I3,'-IH   bJUbNVALUfc    Llr    t4 — HAS   ML!7    B£fcN»  , 

♦  •    OETEPMINDED    AFTER    30    ITERATIONS.     ****♦•/) 

GO    TO    1  ■ 

3    PPINT    2  22,     Ik^ ' 

222  FTRMAT( ////IPX ,•%****  THE  ',13,'  -TH   EIGENVALUE  HAS  NOT  BEEN', 

♦  •  DETEP.MINOFO  AFTEP  3Q  ITERATIONS.  »***»'/) 

GO    TO    1 

A.    PRINT    223 

223  FORMAT (•    ***    MATRIX    P.    IS    NOT    POSITIVE    DEFINITE    DUE    TO', 


"■ «JMPP.JPtP    INPUT    DATA,    r* — =FTjrjTJD   dFF    bukOF.    ***«} 


GO  TO  I 
996-  STOP 
END" 


SO/0218 
SO/0219 


SOZ0220 


37 


CUTINF    507(MfNf*P,nfEPS,wAMT  X  ,"X  r-VHSP  t  »,  ITge  ,  If  LAG) 
ICIT   Ff-AL*3    (a-h,cw) 


TTTTTC 
UPL 
DIMfNSICN    tffMND.ND)  ,0(N'D),X(hD,NOJ  t  ITFS(ND) 


?<V)221 
SQZ0??2 
$070223 


lit i LAi   vm  ;x 

SH7022« 

CALL    SO'TPI  (Nn,N,AB,ntEPStWAMXfX,IFLAG) 

SOZ0225 

CALL    SOZIT(MO,NtAFif  D,EPS,EPSP,  ITEF,  K2fIE6P,WANTX,X,«l  3,  IFLfG) 

S0Z022e 

FfcTUPM 

5UZ0227 

2    FLTURN1 

S0>0228 

3    FfcTURN2 

S070229 

tNU 

5.0/0230 

38 


StWPUTINF    S07T 

implicit  peal*8 
uh;k  ai  waniX — 


0  KNO,  N,  A,  n,  EPSBtWANTXfX,  IFLAG) 
(A-H,0-Z ) 


SQ7D232 

SH70233 
S07O23'' 
SO Z 02 3 5 

SQ7023t 

"STTTTTTT 
SOZ023B 
SQZ0239 

b'J  'U2-i: 

SOZ0241 

sozo2<-2 


DIMENSION    AIND.ND)  f  D(ND)  ,X(  ND,ND) 
3T 


T CIDUCE     A    TU 

IF(N.LE.2)    r,0 
rM2=N-2 

DTI    UU    K  =  l,  K»2 

KP1=K+1 
NK1=N-K-1 
LIU    ibO    L^=l,  r 
L=N-LB 
L1=L*1 

'HA(L1,KJ.  b('.l4 

IF< IFLAG.En.i) 
C    ::::::    SCALING    : 

I  M  a  I  L  »  k  i .  t  i  > .  u  . 


D1AUCNAL,    KbfcP    U   U1ACUNAL. 

n   uo 


KT" 


.  j   i;u    'i!   itu 

00    TO    <83 


SUZ02AJ 
SO? 024* 

SO?n245 


S070247 

SO/0248 


DITJ    GTJ     I  LI    cdi 

1) 

1  ,  K  )  /  A  (  L  t  K  )  ) 
c31  ,iSBi  ,c»J2 
.DO*TSCALE)/(TSCALE-TSCALE*TSC\LE>) 
LE( A,D,LflltN    t    K    , NO.TSC AL E  ,N , w ANTX , X ) 


"WZTZTT 
5^70250 
SOZ0251 
S'l/ J2i2 
S0Z0253 
SO ?D2?4 

tttttzst- 

SQZ0256 

S070257 

i.n/'525B 

S0ZC259 

Sw£2£fcO 

in/.)2oF 

SOZ0262 

S0Z0263 

sons*-. 

SOZ0265 
S0ZQ266 

son?i/ 

S0Z0268 
S0Z02S9 
S07O270 
S0Z0271 
S070272 
SO Z 02 73 


FSf  ALF=D(l.  )/0(L 
TSC  *  Lb=CAQC  U(  L 

lr  I  ISLALTM.in) 

681    IF(RSCALF.GT.(i 
+                    CALL    SCA 
TW3 


TTU-TT 

682  IFIRSCALE.LT. (7 
♦  CALL     SCA 

683  UN!  INUb 

P=DSORT( A(L,K)*- 
F="DSTr,N(P,M  L,* 

C=AIL,K  )/V 

S=ML1,K)/P 
C1  =  C  

5T_3 

IFI  IFLAG.EO.l) 
FDN=C*S*(D(L)-D 


SCALE-l.DO)/ll.DOtTSCALEvTSCALE>) 

LE( A,n,L tLlfN     ,     K    , ND , TS C ALE , N, W ANTX , X ) 


A(L  ,f)«-A(Ll,K)«A(Ll,K)) 
)  ) 


GO    TO    13 
(LI)  ) 


i  =  im  n 

D(L1 )=C*C*T*S*S 
C(L)=D(L  )*T/P( L 


TFJOAB^ED")  .Lb 
P=-EDN/0IL1 ) 
C1=C+S»P 


*D(L) 
1  ) 

.kHSP)    GC   "0    15 


s"i=s-r>p 

15    CONTINUE 

DO    103    J=K,L 


T^  A  { L  1  ,  J  ) 

AlLlf J)=S*A(Lf JJ~C*T 
103    A(L, J)=rl*A(L,J)+Sl«T 


SOZ027- 

"Sl^TTTT 
SOZ0277 


U  =  C1-T+S1'A  (Lit 
A  (LI  ,1.  1)  =S*T-C  + 
Mlf  L)=Cl*A<L,L 


LI) 

ML1  ,11  ) 
(♦S1MJ 


04    A(J,L  l)  =  S*T-C  A(  J,L  1  ) 

I  F(  .NQT.WANTX)     f-C    TC     10' 
DP    10    J=l ,N 


T  =  X( J,L  ) 

X( J,L )  =  C1»T  +  S1  vxrj,Li) 
10    X( J,  I  1) =S*T-C*X( J,Ll) 


S;0  7O?79 
Sr,Z02S0 
S07  T291 


DO     1C4    J  =  Ll,r; 

SOZ  )162 

T=A( J,L ) 

SO/0233 

MJ,L)  =  C1»T  +  S1*MJ 

,L1) 

S07  128- 

SO  7 12  ii 
SO  702 36 
< 070  2d  7 


SQ7  )2Jb 
S0702J9 
SO  7 0200 


105    CCN7IMIE 


prnzn 


39 


~~T5U" CTRTTNUr 

160    CONTINUE 
170    CPNHNUb 


TTTTJR^T 
END 


SQ70292 
SOZ0293 
SQ7  02K 

S0Z0296 


Uo 


c 

SUBROUTINE    SCV!T(Nn,N,A,U,EPS,  E«>SB,  ITE*,*, 
IMPLICIT    REALMS     (A-H,C-Z) 

IERR,WANTX,X,*,IFL*G) 

SQ7>>297 
S0Z0298 
SOT0299 

ICG ICAL    WAMTX 

DIMENSION    A (NO, NO) t D( NO ) , X( ND, NO) 

DIMENSION    I  TFP  (NO 

SOZ0301 
SO70302 

c 
c 

INITIALIZE     ITFF,    COMPUTE    EPSB 
LS  =  1 

bCJZ0J03 
SOZ0304 
SOZ03">5 

"    V 

BNORM    =    0. 
DC    185    1=1, N 

b'JZWJi. 
SOZ0307 
SQZ0300 

I  I t«  (1  )     =    0 
BNI=DABS(D( I  )) 

IP     (BM.GT.RNORM)    FNOPM    =    PNI 

sozjj')^ 

S0ZD310 
S0Z0311 

18b 

U'\i  Il^'Ut 

IF     (PNCRW.EO.O.)    BNCPM    =    EPS 

EPS6=EPS 

bOZ'JJU 

SOZ0313 
S0Z03K 

c 

REDUCE    A    TO    DIAGONAL,    KEEP    B    DIAGONAL: 

M    =    N 

% 

S'1703I! 
S0Z0316 
SQZ0317 

c 
c 

700 

lh      IM.Lb.LSJGU     \U     390 

CHECK    FTP   CONVERGENCE   OR    REDUCIBILITY 

by; jjid 

S0ZJ319 

S3Z0320 

c 

DO    220    LP=1,W 
L    =    M+l-LB 

SOZ0322 
SOZ0323 

c 

LMl  =  L-l 

IF    (L.EO.LSJOO    TC    260                      --- -    ... 

50 Z J  324 
SQZ0325 
S0Z0326 

IF{DARSCA(L,LMl)  J.Lfc.fcPSJ    GLi    TC   230 
TAU1=A( LM1.LM1 ) 
TAU2=A(L ,L) 

SQZ0328 

SOZ0329 

GiMMA=i JL.LM1)                                                       '     ' 
IF(  IFLAG.EO.l)     GO    TO    219 
DLDL1=D( L)*D(LM1 ) 

SrtZ}33fl 
S0Z0331 
SQ70332 

IHDLDL1.L1  .1.)     ktlUANl 
Dn=DSORT(DLOLl) 
IFtDLDLl.GT.0. 1    GO    TC    2H 

SP73333 

S0Z033^ 
SOZ0335 

214 

GC    TC    220 

CONTINUE 

LP1=L*1 

SOZ0337 
SOZ0338 

DESP1=16.D-1? 
IF(DD.GT.DESP1 )    GO    TO    218 
DESP2=16.0+03 

?^Z0i39 
SQZ0340 
S0Z")3'  1 

DESP3=0£SP2-DESP2 

DO    217    1=1,2  .  S07T343 

LP1VI  =  LP1-I S0Z3344 

DILP1«I  )=0(LP1WI)'DESP3  SO703AT 

IF(LP1MI.GT.LSJ     AiLPl^I  ,  LPU'I-1)  =  A(LP1MI,LPIMI-1)*DESP2  SOZ034*- 

A(l  PI  Ml  ,LP1M!  IsAtLPlMItLPl*!  )»0ESP3  SOZCH',7 


216 

IK  .f.OT.VANTX)     GO    TO    217 

DO    216    J=1,N 

X.(  J,lPlVI)=X(JtLPlMI)*DESP2 

SO 1 0348 
SOZ33^9 

SO703=O 

217 

CONTINUE 

ML.L^l  )  =  A(l .,L«1I*DESP2 

IF(LPl.LF.N)    A«LP1,L)=A(LP1,L) <0cSP2 

SO/0  351 
SQZ1352 
S"70353 

218 

c 

C 

CCNTINUE 

FIFST    STPPINr.    rriTFRIrN    : 

SO/035^, 
S1Z03^f 
$0/T33<L 

TAU1  =  M  LM.L^'l  l/OlL'ni 

SCZ0357 

Ul 


fJU2=A( LtU/niLI 

GAM^A=A(ULM1)/DD 

219    CONTINUE 

Z T  =UArib(  !AUimUrfS(  IAU2  J 

EPSAA=EPS  *0T 
JF<DABS(GAMMA)  .LE.EPSAA) 


GO  TC  230 


5070356 
SOZ3359 
SQ703iO 
SO  Z  ^  -3  c  1 
SQ7D362 
SQ70353 


c 

i0/03o4 

c 

SECOND    STOPPING    CRITERION    :. 

SQZ0365 

CC=0APS{TAU1*TAU2) 

SQZ0366 

hPSI=bPb*+A 

5rV  J36/ 

IFIGG.LE.EPSM    GO    TO    220 

S0ZO368 

TAU12  =  <",AP"A»t2/DAHS (TAU1-T AM 2 ) 

S0703o9 

IH  IAU12.Lt.tHi     J     GU     IU     III 

SOZJ370 

GC    TO    220 

SQ71371 

777 

GS=GAMMA/GT 

SOZ0372 

LPSM    rfcPS    *(1*.U0J 

SCZJ3f3 

IF(GS.LE.EPSM    1    GO    TO    230 

S0ZO37<- 

220 

CONTINUE 

SOZ03  75 

230 

£IL,L-1)    =    0. 

S0Z0376 

IF    <  L.LT.M       )    GO    TO    260 

S020377 

M    =    L-l 

SOZ0378 

rc  70  200 

!>yZ0if9 

c 

S0Z0380 

c 

CHECK    FOP    S^ALL     TOP    OF     R 

SQZ0361 

c 

S0Z03J2 

260 

IF     <DARS<    0(L)     I.GT.EPSB)    GO   TO   300 

SOZ0333 

D(U  =  0.000 

S0Z0384 

LP1=L*1 

?a?.")jd.e 

LP2=L+2 

SOZ0366 

lP(LPl.GT.M)     LP1=M 

S0Z0387 

1HLP2.GT.M)    LP2  =  M 

5QZ0338 

ll=LPl 

SQZ0389 

IF(A(LP1,L  )  .F^1.  J.OOO)    r,Q    TO    280 

S0ZO390 

P=-A(LPl,L)/A(L,L) 

S0Z0391 

IF(DARS(P) .LE.l.OOJ     GC    TO    27C 

S0Z0392 

IMsiniNTfDLrGIDABSIP) )/DLOG(  16.D0J )  +  l 

S0ZO393 

6*16. bo**  im 

SQZ0  39'' 

A<L,L)=Q*Q*ML,L) 

SOZ0395 

A  (LPl,L>  *0*  f  (1*1, L) 

SnZ0396 

IF( .NOT.WANTXI    GO    TO    269 

SOZ0397 

DO    265    1=1, N 

S0Z1398 

265 

X(I,L)=0*X(I,LI 

S0Z0399 

269 

CCNTINUE 

S0Z04)0 

P=-t (LPl,L)/A(Lt I  ) 

SOZO-01 

270 

CONTINUE 

S0Z0vO2 

A( IP1,LP1)-A (LP1.LP1 )*P-A(LPl,LJ 

SOZO-'-03 

IF( .NPT.WAMTX)     GO    TO    279 

SQZ040* 

DC    27=;     1  =  1  ,N 

SQZ0AQ5 

275 

X(  I  ,LP1)  =  X(  I,LP1 ) ♦ P *  K { I ,L) 

SOZO-»Ofc 

279 

CCNTINUE 

SOZ0*t07 

280 

I  =  L  P  1 

S^Z  •■)••  c  8 

LS=LS*i 

S^ZJ09 

GC    TO    2  30 

SQZOUO 

c 

SOZ0«Ul 

c 

PEGIM    ONE    SOZ    STFP,     ITERATION    STRATEGY 

SOZ0M2 

c 

SHZO  >13 

300 

V  1  =  M-  I 

SQ?OU<- 

LP1=L*1 

$OZV  15 

LP2=L*2 

SOZOVlc 

IF(LPl.GT.M)     LPl=v 

S0Z3M7 

IKLP2.r,T.M|  LP2  =  '" 


S0Z0-13 


k2 


ITER<w)=ITER(M)*l 
IF(  ITER  <»').F  0.1)    GO    TC 
Ml)}     .LI. 


'TFTn*R5(AIM, 
1F(ITER<M)     .LE 

IEPf-=M 
RtTUPMl 


305 


30)      r,o 


TTnO*CLDU 
TO    305 


CJ   73   305 


SrV)'*20 
SO/0^21 


SGZ0'.22 
S0Z0t23 
S0?0V2^ 


S0Z042: 
SOZ0^26 
SOZO-'27 


305    CONTINUE 


S0Z0«»2(J 
S070«*29 
SQZO'OO 


TO    200 


TUTFTT 
S0Z0'+32 
SOZ0<>33 


2»czut3- 

S0ZOA35 
SOZO- 3e 


INV) )    *DARS<  A(M,M1)) 


"suzn-3/ 

S0Z0v38 
SOZ0<-39 


IFJP     .LT.     0.00)     R=-P 
SHIFT=G2+EL/(P+R) 
::::::    END    OF    FINDING    S H I*- T 


'suznuo 

SO?OM1 

S0Z0^42 


IF( IFLAG.EO.l) 


TT^ ML.L) 

X2  = 
IF(X2.EO.O.O0) 


GO    TO    583 
-    JjHIKMilL) 
A(LPltL) 
GO    TO    2  80 


S0Z0+4* 
SQZ0t«r5 


VJ'.O'Tt 

S0Z0^7 
S0  7>,-o 


n 


33 


::     SCALING    : 

IF(  Xl.FO.0.00)    r,c 

TSCiLFsClLJ/lKLi'l  1 

TSCALE=DAPS (X2/XI) 

IFtTSCALE-1  .DO)    =81,583,582      

I  1.  Jd+'UALU/l  I  .jcALc-i^ALb-IWltlJ 


SQZO^O 
SQZ0451 


58L     IFiL'SLALF.?;". 

+                    CALL    SCALE(  A,D,L,LP1  ,  LF2  ,L,  ND,  TSCAL  E,  N,  WAN'TX  ,  X  ) 
GO    TO    583 
5  82    IFlKiCALL.L'.l  I  SL4L  l-1  .  DO  )  /  (  1 .  L'O+TVJ  tLt  -i  "C  ALc  )  J 

♦  CALL     SCALE(A,D,L ,LP1 , LF2 , L , NO t TS C ALE ♦ N, WANTX , X ) 

5  83    CONTINUE 


SUZrt-52 

sqzo*53 
son45* 


SQZ3"56 
Sn?n-T57 


TrT>T^TnrrrrTT(MtM-i)) 

Xl=  A(L,L)    -    SHIFT+DIL) 

X2  =  MLP1  ,L) 


SOZViB 
S0Z0-'f9 

sozrvoo 


SQZ  i-bl 
S0Z0462 
S0Z0'-c3 


PFP 
CAL 
IF( 
DC 
KM1 
Kl  = 
K2  = 
IF( 
"I  f  ( 
IM 
IK 


1F( 
PSr 


681 


TfC 
IF( 

If  ( 


FTP 
L  T 
I.  PI 
3  60 
=  H  — 

K+r 

K+2 

K^l 

kTT 

.»  IK 
I  Fl 
SrA 
A(K 
/■IF 
ALE 
T5C 
P  Sf 


M    THE    SPECIAL    TRANSFORMATION    PRODUCED    FR0M    THE    ORIGIN    SHIFT: 
ftArSF<A,L),L,L,LPl  ,  I  P2tU/N:TX,X,Xl,  X2 , N, *1 , M,ND ,  EPS R, I  FLAG) 

.gf.v)   c,:, 

K=LP1,M1 

1 


TO    3  30 


S0ZDvO<r 
SQZ0-65 
S0ZO"6c 


•LT.L)  K"l=L 


S0Z0468 
SOZ0v5  9 


OT.^)    K2  =  v 

l,*Ml )     ,E0.    O.CO)     CO    TO 

AG.EO.l)    r-7    TO    683 


360 


SO/0*  ZO 
SQZO-:  71 
Sn?DW2 


L  I  KG 

,K^1) 


.EC.n.OI)     GO    TO    683 

/nun 


S?Z0~Z3 

sozn^ 

S0Z-)v7  = 


=  D  Ah  S  (  /  (  k  1  ,  «'*  1  )  /  'l  (  u  f  k  '■'  1 )  ) 

flLr-l.TP)     -81,683,632 

ALF.0T.(l.:K)f7SCALfc  )  /  (TSCALE-TSC  *  LF.»TSC  ALE  H 


Snz0v76 
SOZC-- 77 
S0'.v»78 


S"70--7» 


CALL     SC/LF  (  A  ,0,K  ,M,K2,  K'-'l,NiJ,"bCuL^f  N.wANTX,  <) 


k3 


GO   TO  6  81 

682     IF(PSCALE.LT.(TSCALE-l.DO)/(l.CO«-TSCALF*TSCALE)  ) 

♦  CALL     SC/LF(  A  .D.KtKl.KZ.K^l.ND.TSCALE.N.WANTX.X) 


SQ7O430 
SQ^OhSI 
SQZO-,8  2 


c 

c 

683 

CUN1  INLIb 

PERFOCV    TRANSFORATIONS    TO     ZERO    A(K 

♦  1,K 

-I) 

C    A(K- 

■It 

K  +  l) 

AND 

S^'o' 83 
S0Z0A84 

S070-65 

L 

KELP"    H    DIirJ'N*L     : 
X1=A(K,KM1) 
X2=MK  1,KM1  ) 

b0/j--3d 
50*0487 

SOZO<,38 

360 
330 

CALL    TEA^iSr  (  A,D,KM1  ,K,nf  K2,WJKTX,X 
CONTINUE 

CONTINUE 

.111 

X^,N 

iMl,H 

H\) 

»bP59 

iIFLASI 

SC/0-89 
SQ 70490 
SOZ0A91 

C 

C 

END   FAIN   LLCP 
GO    TO    200 

SQZ">493 
SQZ049' 

c 
c 
c 

END    ONE    SQZ    STEP 

SQZ0496 

sozo;-97 

390 
555 

(TNI  INUfc 

WP.ITE(6,*?5?  )          UTERd!  )  ,II«1,N) 

FOPMATM    *»■*+**    ITEMI  )  =  •/,  (2015  )  ) 

SOZ0<t99 
SO  7 05 30 

FETURN 

SO  70;  Tl 

END 


S0Z0502 


kk 


c 

SQ7O503 

SUBROUTINE       TR A\'SF( A , D, KM1, K ,K 1 ,K2,WANT X, X, XI , X2. N,  M I ,M ,ND, EP5B, 

S0Z0«534 

f                                                  IFIAO) 

SQ 7^505 

IMPLICIT    P.E£l>8    (A-H.C-ZJ 

50Z05Ce 

LOGICAL    WANTX 

SOZ0S37 

DIMENSION  AtNOtNDN  0(N01  tX(NDt.NO) 

SO'03C8 

c 

50Zv350«J 

c 

FIND    AND    PERFORM    TP  ANSFOCMATITNS    TO    ZERO    THE    UNWANTEO   NONZERO 

S<V0510 

c 

ELEMENTS    IN    MATRICES    A    AND    R     : 

SO70511 

I 

i>UZ  J- 12 

P  =  DSf>RT(xi*XH-X2*X2  ) 

S070513 

R=DSIGN(R,X1) 

SO  70 5 14 

l=Xl/K 

SUZJal^ 

S=X2/° 

SOZO'lt 

CS=C*S 

S073317 

LL  =  L*«. 

5073318 

SS=S*S 

SOZ0519 

Cl  =  C 

SO  70  520 

S>1  =  5 

5CZ1i2l 

IF(IFLAG.EO.l)     GO    TO    7 

S0ZO522 

EDN  =  rS"*  (D(K  l-D(Kl)  ) 

S070523 

T=D«K1) 

SOZfliZA 

D(Kl)=CC*T+SS+D(K) 

SOZ0525 

D(K)=D(K)*T/0(K1 ) 

SQ70526 

IFJDJBSIEDNJ.LE.E^PJ.  C^   TCI    7 

5070^27 

P=-EDN/D(K1 ) 

S070526 

C1=C*S*P 

SO 70-29 

S1=S-C*» 

S0M330 

7 

CONTINUE 

S070531 

DC    21    J=KV1,K 

SO70332 

T=MK1,JJ 

SQJ0333 

A(Kl  ,JJ=S*A(K, J)-C  +  T 

S07053^ 

21 

* (K, J)=C1«A(K, J)+Sl*T 

SO  ?0«?3E 

IFU.G7  . k M l  J    Ain,^ll=0. 

507053c 

U=C1*T+S1*A(K1,K1) 

SQ703-37 

A(K1  ,Kl  )  =  S*T-C- A(K1,K1) 

S07O538 

A(K,K)=CI-A(K,<)+S1'U 

5Q70  5jS 

T=A«K1,K1) 

SQZ0540 

A(K1  ,K1  )=S*MK1,K)-C*T 

SOZO^l 

A(Kl  ,K  )=ci*  MKt,t<l+SITT 

SQ703''2 

IF(K2.E0.K1 )    CO    to    23 

S070543 

A(K2,K)=Sl*AfK2,Kl) 

S0Z0544 

A(K2  ,K1  )=-C*A<K2,Kl) 

SQZ05-5 

23 

CONTINUE 

SQZ05*6 

IF( .N0T.WANTX1     GO    TG    30 

S0705',7 

CO    25    J=1,N 

SQZ03'>8 

T=X( J,K1) 

SQ705-VS 

X(  J,*l )=S' X(J,K)-r»T 

S370^CP 

X(J,K |=C1*X(J,K)+S1*T 

SOZO?i 1 

23 

CCNT INUE 

S070332 

30 

CCNTINUE 

SO  7 0  55 3 

RETU 
END 


s  rj 


SQ70  554 

SOZ05?5 


k5 


5 

SUFFCUTINE    SCA LP ( t , 0, K , Kl ,K2 , K«l , ND.TSCAL E ,N ,WANTX , X ) 
ILLICIT    P.EAL-»8(A-H,0-Z) 

SOZOC57 
SO  >  0^8 

c 

LT&ical  utwrv. 

DIMENSION    A(NOtMO) »D(NO) ,X(NC»ND) 

so7oc:9 

SOZOC60 
S0705M 

c 

c 
c 

THIS    SU3FGUTINE    F  T  N  n  5   THE!   SCALIKG  FACTOR    FCR   STA6ILIHNC   THE" 
ELEMENTARY    TRANSFORMATION    : 

SQZO'62 
S0Z0363 
S(3?05£' 

THl<",H=l4,.DO-le.no 
TL0W=1. DO/THIGH 
IF(TSCALE.LT.TLOW)     GO   T0    888 

SOZ0566 
SCZ0*67 

888 

IF(TSCALE.LE.TH!'-,HJ     GC    TO    889 
SCALE1=DS0FT(D(K1) ) 
?CALb2=0SOPT(D(K) ) 

SQZO=f 9 
SO7.0Z7C 

r 

:::    NOW    iJDniTF   THE    AOWS    AND   CCLUWM5    K.K1    CF    A    AND   D    : 
0(K)  =  D(K  )*0(K1 ) 

p(ki  )=n(K) 

5SZ0?71 

SOZC572 
SOZO'73 

998 

f  CN7INUE 

DO    12    J»KM1,K 

MK,  J)*$C*LE1*A(K,J) 

SQZ0575 
SOTOf 76 

12 

AlKl  ,  JI-SCALF.2'  MK1,  J) 

A{K  ,K)=SCALE1*A<K,K) 

MKl,Kl)«SCJLE2*A(KlfKl)vSCALE2 

SQ70;77 
SQ70373 

SOZO'79 

A(KlfK)»SCALEl*AiKl,"1?l 
IF (K2.E0.K1 )     GO    TO    13 
MK2»K1)*SC*LE2*A(K2.K1) 

so'.osan 

SQZ0381 
SQZ0i82 

13 

CONTINUE 

IF< .NCT.WANTX)    GP    TO    20 

00    1*    J*1,N 

SQZ0583 
SO  Z  05  8'- 

S0Z05dc 

1* 
20 

X<J,K)  =  SC.*LE1«-X(J,K) 
X( J,K 1)=SCALE2<  V( J, Kl) 
CCNTIK'UE 

SOZ0566 
S0Z0587 
S0ZC588 

889 

GC    TO    1818 

SCALEl=l.DO 

SfALE2*l.0n/TSrALE 

S0Z0:o9 
SOZ0590 
S070391 

D(Kl)*ISCALE2T«;CALE2K0(Kil 

S0^0=i92 

rr    TO    998  SOZJ593 

1618    CCNTINUE . S0Z0C9* 

RETURN  SQZ0s9f 

END  SQZ0596 


U6 


SUPPHITINE    CODER  INiNDfD»X) 
IMPLICIT    RFAL»8( A-H,CW) 
"DIMENSION    x"l,\iD,MJ)  ,!1IND] 


SO?0C97 
S070-98 
SO?0r99 

sanooo 

S07G601 

sojrwv 

5^/3=03 
SOZPcOf 

"STrrj^jr 
SQZ06C7 
S0ZC6GB 

SQZ.OclO 

S07H613 
SQZOal* 

SQ/Otl6 

SO' 1617 

SOZ^lfi 
SQ/0619 
SQZO620 

tttottt 


THIS    SlBf-miTINE    CR-OEPS    THE    EIGENVALUES    OF    MATRIX    B     IN    ASCENDING 
TPDEP    S3   that    THE    4ETTE;.   RESULTS    I U    THE    SECOND    CiU    OF    Wl   wlLL 
BE    GUARANTEED. 


DC     20    K=l,f, 
DMIN=1.  00+-38 


JfIN=K  

T5TJ    10    J  =  K,N 

IF(DMIN.LE.D( J) )    GO    TP    10 
DMIN=D( J) 


DMIN=D( J) 

JWTTTT] 

10    CONTINUE 


T  =  D  (  J  M  I  N  )  

TJUMTNi     »=DJKJ 

D(*)=T 

OP    2  0    I  =  1 , M 

T  =  X(I,JM!M 

X(I,  JMI'j)=X(I  ,K) 


20    X(  I,K)  =  T 

PFTUPN 

END 


SOZOo22 


hi 


"STiFROUTIfE    X7rX(ND,N,C,X,D,Wl) 
IMPLICIT    RE*L*8U-H,0-7I 
DIMENSION    C(ND,ND)  ,X(ND,*'D),Wl(NO),0(NO) 


S070&23 
son  --24 
SQZPc2c 


THI  S 
THIS 


SIJPFOUTINE 
SUROf.uTINE 


UPDATES 
IS  ALSO 


THE  MATRIX  A 
USED  TO  PEP. PC 


AFTER  FIRST  CALL  OF  SOI. 

RM  ANY  MATRIX  MULTIPLICATION 


S'3  7  jc  2£ 
SOZ0627 


■urTHrrrp — XT   *   Z  *  X. 


S0/0c29 
SQZ0630 
SQ70631 


c 

c  —  c  *  x 


Hi i  yo  j  =  i,\ — 

JM1=J-1 
IF(J.LE.l)    GO 


5TJ7TToT2 
S070t33 

SO?1t34 


tO    * 


DP  3" 
KM1  =  K 
KP1=K 


K=l,JVl 

-1 

+  1 


TTTTTTr 
S070636 

<OZ0637 


Ullk)  =  D(K)     y (K  f J  ) 

DC    1    I  =  KPL,r.' 
1    Wl  (K  )=W1  (K)  Kf  ,1  )«X<  I,  J) 


"snTTTIB" 
SQ70h39 
SOZ0640 


TMK.LE.l)    r,.j   TC    3 

DC    2     I  =  1»KM1 
2    W1(K )=W1 (K)+r( I ,KJ+X(I , J) 


nJTJFTI 

S070r42 
SnZ0ti3 


3  COT! 

4  CCNTI 
DC    30 


RTJE 

NUE 
K  =  J,N 


SQ?  064-5 

Sn?06'-6 


KPl  =  K 
Wl  (K) 


=1 

=D(K) *X(K, J ) 


ToTTtTT 
SO Z 0*48 
S070*-9 


IPJK.tO.M     r-C    TO    I  if 

DC    10    I=KP1 ,N 

V1(K  )-Wl  <K)+C(K,  I  )«-X(I  ,  J) 


S07O6=i 
SOZ  )652 


10 


15    CCNTI 
IF(K. 

DO    20 


*TJZ 

EQ.l)    OC    to    30 
I=liKMl 


T<T7DTTT 
SQZ0i5* 

SOjTQc5? 


20    wllKl 

30   CCNTI 
C XT    * 


=  wUkJ+C(!,KJ«xII,JJ 

nue 

(C    *    X )    : 


S0Z0657 
SOZOi5:8 


DC    56 
C  (  I  »  J 

Dr_L2 
TTlt  J 

ccnti 

FETUP 

END 


I=J,N 
1=0. 
K  =  1 ,  N 


S^ZO-59 
SQZ0*60 

snzo*>6i 


50 
90 


)=C<  I, J  l*X{K, I  )*W1(K) 

NUE 
N 


SQZ0*c2 
S0ZC663 
Snzof6<- 
SQZ0665 


**c. 


hi 


SUPRCUTUiC    Xf("X<ND,N,C,X,P,WU 

IMPLICIT    RE£L*8( A-H,0-7) 

DIMENSION    C(ND,ND)  ,  X  (ND,*'D)  ,  W  1  (NO  ) ,  Q(ND  ) 


S070&23 
SQ'0624 
SQZPc2c 


53/ jc2£ 
SOZ0627 

SQZn*2e 


THIS    SIJPFniJTINE    UPDATES    THE    MATRIX    A    AFTER     FIRST    CALL    OF    SO?.. 
THIS    SUBROUTINE    IS    ALSO    USED    TO    PEP.FQR.M    ANY    MATRIX    MULTIPLICATION 
OF    THE    FCTy XT 


TT 


5070r2<5 
S0ZP630 

sorot3i 


~s~T 


C    *    X    : 


SUMeJ2 
S070633 

S070S3^r 


li'j  yu  j  =  i,n 

JM1=J-1 

IF(J.LE.L)    OC    TO    * 


TTTTT3T" 
SQZ0636 
SOZ0637 


00   3 

KM1  =  K 
KP1=K 


K=1,JV1 

-1 

+  1 


U1(KI=D(KJ     X|k,JJ 

DC    I     I»KPl,N 

W1(K)=W1  (KltCIM  )«X(  I, J) 


SQ70639 
SnZOc-':0 


TF(K.LE.l)    f,.j    il    3 

OC    2     I=1,KM1 

Wl(K)=WKK)*C(ItK)*X(I,J) 


S070642 
SnZ0ti3 


CCNTI 
CCNTI 
DC    30 


RUE 
NUE 
K  =  J,N 


SQ?06^E 

Soroe'-ft 


KH1  =  K 
KP1  =  K 
W1(K) 


=1 

♦  1 

=D(K) *X(K,J ) 


S070t47 
SQZ0<-48 
SOZO^-9 


IF(K.tO.N  J    r.i,    TO    lb 

DC    10    I=KP1 tN 

V'K*  )=W1  (K)+C(Kf  I  )*X(I  ,  J) 


TOTTTTO- 

S0ZO6;l 

SOZ  )652 


10 


CCNTIMJE 

IF( K.EQ.l)    OC    TO    30 

DO    20    1=1, KVi 


"PT70TTT 

sozots* 

S^ 70c 55 


20     W1(KJ=WI{K)*C"(!,K)«XII,J) 

30   CCJ^'TINUE 
C    XT    +    (C    *    X )    : 


ToTnTTT 

S0Z0657 
SQZ0*58 


I=J,N 
)=0. 
K  =  1,N 


S070i59 

sozr*60 
sozo*>6i 


DC  50 
C(  I,J 
DO    5  0 


TTi,  J 

CCNTI 
FETUP 
END 


)=C(  I, J  >  +  X(K,  I  )*W1(K) 

NUE 
N 


50 
90 


S0Z0*c2 
S0ZC663 
Sfl?0*><s<- 
S0Z06c5 
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